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Abstract: We determine the solar neutrino fluxes from a global analysis of the solar 
and terrestrial neutrino data in the framework of three-neutrino mixing. Using a Bayesian 
approach we reconstruct the posterior probability distribution function for the eight nor- 
malization parameters of the solar neutrino fluxes plus the relevant masses and mixing, 
with and without imposing the luminosity constraint. This is done by means of a Markov 
Chain Monte Carlo employing the Metropolis-Hastings algorithm. We also describe how 
these results can be applied to test the predictions of the Standard Solar Models. Our 
results show that, at present, both models with low and high metallicity can describe the 
data with good statistical agreement. 
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1. Introduction 

The idea that the Sun generates power through nuclear fusion in its core was first suggested 
in 1919 by Sir Arthur Eddington who pointed out that the nuclear energy stored in the 
Sun could explain the apparent age of the Solar System. In 1939, Hans Bethe described 
in an epochal paper [1] two nuclear fusion mechanisms by which main sequence stars like 
the Sun could produce the energy necessary to power their observed luminosities. The 
two mechanisms have become known as the pp-chain and the CNO-cycle [2]. For both 
chains the basic energy source is the burning of four protons to form an alpha particle, 
two positrons, and two neutrinos. In the pp-chain, fusion reactions among elements lighter 
than A = 8 produce a characteristic set of neutrino fluxes, whose spectral energy shapes 
are known but whose fluxes must be calculated with a detailed solar model. In the CNO- 
cycle the abundance of 12 C plus 13 N acts as a catalyst, while the 13 N and 15 beta decays 
provide the primary source of neutrinos. 

In order to precisely determine the rates of the different reactions in the two chains, 
which are responsible for the final neutrino fluxes and their energy spectrum, a detailed 
knowledge of the Sun and its evolution is needed. Standard Solar Models (SSM's) [3-9] 
describe the properties of the Sun and its evolution after entering the main sequence. The 
models are based on a set of observational parameters (the present surface abundances of 
heavy elements and surface luminosity of the Sun, as well as its age, radius and mass) and on 
several basic assumptions: spherical symmetry, hydrostatic and thermal equilibrium, and 
equation of state. Over the past five decades the solar models were steadily refined as the 
result of increased observational and experimental information about the input parameters 
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(such as nuclear reaction rates and the surface abundances of different elements), more 
accurate calculations of constituent quantities (such as radiative opacity and equation of 
state), the inclusion of new physical effects (such as element diffusion) and the development 
of faster computers and more precise stellar evolution codes. 

Despite the progress of the theory, only neutrinos, with their extremely small interac- 
tion cross sections, can enable us to see into the interior of a star and thus verify directly 
our understanding of the Sun [10]. Indeed from the earliest days of solar neutrino re- 
search this test has been a primary goal of solar neutrino experiments, but for many years 
the task was made difficult by the increasing discrepancy between the predictions of the 
SSM's and the solar neutrino observations. This so-called "solar neutrino problem" [11,12] 
was finally solved by the modification of the Standard Model of Particle Physics with the 
inclusion of neutrino masses and mixing. In this new framework leptonic flavors are no 
longer symmetries of Nature, and neutrinos can change their flavor from the production 
point in the Sun to their detection on the Earth. This flavor transition probability is 
energy dependent [13-16], which explains the apparent disagreement among experiments 
with different energy windows. This mechanism is known as the LMA-MSW solution to 
the solar neutrino problem, and affects both the overall number of events in solar neutrino 
experiments and the relative contribution expected from the different components of the 
solar neutrino spectrum. Due to these complications, at first it was necessary to assume 
the SSM predictions for all the solar neutrino fluxes and their uncertainties in order to 
extract reasonably constrained values for neutrino masses and mixing. The upcoming of 
the real-time experiments Super-Kamiokande and SNO and the independent determination 
of the flavor oscillation probability using reactor antineutrinos at KamLAND opened up 
the possibility of extracting the solar neutrino fluxes and their uncertainties directly from 
the data [9, 17-22]. Nevertheless, in these works some set of simplifying assumptions had 
to be imposed in order to reduce the number of free parameters to be determined. 

In parallel to the increased precision of the SSM-independent determination of the 
neutrino flavor parameters, a new puzzle has emerged in the consistency of SSM's [23]. Till 
recently SSM's have had notable successes in predicting other observations. In particular, 
quantities measured by helioseismology such as the radial distributions of sound speed and 
density [5-8] showed good agreement with the predictions of the SSM calculations and 
provided accurate information on the solar interior. A key element to this agreement is the 
input value of the abundances of heavy elements on the surface of the Sun [24] . However, 
recent determinations of these abundances point towards substantially lower values than 
previously expected [25,26]. A SSM which incorporates such lower metallicities fails at 
explaining the helioseismological observations [23], and changes in the Sun modeling (in 
particular of the less known convective zone) are not able to account for this discrepancy [27, 
28]. 

So far there has not been a successful solution of this puzzle. Thus the situation is 
that, at present, there is no fully consistent SSM. This led to the construction of two 
different sets of SSM's, one (labeled "GS") based on the older solar abundances [24] imply- 
ing high metallicity, and one (labeled "AGS") assuming lower metallicity as inferred from 
more recent determinations of the solar abundances [25,26]. In Ref. [9] the solar fluxes 
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corresponding to such two models were detailed, based on updated versions of the solar 
model calculations presented in Ref. [8]. These fluxes were denoted as "BPS08(GS)" and 
"BPS08(AGS)", respectively. In a very recent work [29] an update of the BPS08(AGS) 
solar model has been constructed using the latest determination of the compositions [26] as 
well as some improvement in the equation of state. For what concerns the overall normal- 
ization of solar neutrino fluxes, the predictions of this new model are very close to those of 
BPS08(AGS). 

In this work we perform a solar model independent analysis of the solar and terrestrial 
neutrino data in the framework of three-neutrino masses and mixing. The aim of this anal- 
ysis is to simultaneously determine the flavor parameters and all the solar neutrino fluxes 
with a minimum set of theoretical priors. In Sec. 2 we present the method employed, the 
data included in the analysis and the physical assumptions used in this study. The results 
of the analysis are given in Sec. 3, where we show the reconstructed posterior probability 
distribution function for the eight normalization parameters of the solar neutrino fluxes. 
We discuss in detail the effect of the luminosity constraint [30] as well as the role of the 
Borexino experiment and its potential for improvement. In addition, we use the results 
of this analysis to statistically test to what degree the present solar neutrino data can 
discriminate between the two SSM's. Finally in Sec. 4 we summarize our conclusions. 

2. Data analysis 

In the analysis of solar neutrino experiments we include the total rates from the radio- 
chemical experiments Chlorine [31], Gallex/GNO [32] and SAGE [32,33]. For real-time 
experiments in the energy range of 8 B neutrinos we include the 44 data points of the elec- 
tron scattering (ES) Super-Kamiokande phase I (SK-I) energy-zenith spectrum [34], the 34 
data points of the day-night spectrum from SNO-I [35], the separate day and night rates for 
neutral current (NC) and ES events and the day-night energy-spectrum for charge current 
(CC) events from SNO-II (a total of 38 data points) [36], the three rates for CC, ES and 
NC from SNO-III [37], and the 6 points of the high-energy spectrum from the 246 live days 
of Borexino [38] (which we denote as Borexino-HE). 1 Finally, we include the main set of 
the 192 days of Borexino data (denoted as Borexino-LE) [40] in two different forms: in 
one analysis we use the total event rates from 7 Be neutrinos as extracted by the Borexino 
collaboration, while in the other we perform our own fit to the Borexino energy spectrum in 
the region above 365 keV (corresponding to a total of 160 data points). Full details of our 
Borexino data analysis are presented in Appendix A. In the framework of three neutrino 
masses and mixing the expected values for these solar neutrino observables depend on the 
parameters Am^, 6*12, and 6*13 as well as on the normalizations of the eight solar fluxes. 

1 We have not included here the very recent results on the low energy threshold analysis of the combined 
SNO phase I and phase II [39]. These results provide information on the 8 B and hep fluxes and show 
no major difference with the results from their previous analysis, hence we expect that they will have no 
important impact on the results of the global analysis here presented. In particular we notice that their 
best fit determination of the 8 B flux as well as of the oscillation parameters are in perfect agreement with 
our results. 
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Table 1: The reference neutrino flux tf>\ used for normalization, the energy ct^ provided to the 
star by nuclear fusion reactions associated with the i th neutrino flux (taken from Ref . [30] ) , and the 
fractional contribution of the i th nuclear reaction to the total solar luminosity. 



Besides solar experiments, we also include the latest results from the long baseline 
reactor experiment KamLAND [41,42], which in the framework of three neutrino mixing 
also yield information on the parameters Am^, #12, and #13. In addition, we include the 
information on #13 obtained after marginalizing over Am^, 623 and 5 CP the results from 
the complete SK-I and SK-II atmospheric neutrino data sets (see the Appendix of Ref. [42] 
for full details on our analysis), the CHOOZ reactor experiment [43], K2K [44], the latest 
MINOS disappearance data corresponding to an exposure of 3.4 x 10 20 p.o.t. [45], and 
the first MINOS — > v e appearance data presented in Ref. [46]. Details of the oscillation 
analysis of these observables will be presented elsewhere [47]. 

In what follows, for convenience, we will use as normalization parameters for the solar 
fluxes the reduced quantities: 

i 

with i = pp, 7 Be, pep, 13 N, 15 0, 17 F, 8 B, and hep. The numerical values of are 
conventionally set to the predictions of the BPS08(GS) solar model and are listed in Table 1. 
With this, the theoretical predictions for the relevant observables (after marginalizing over 
Am| 3 , 623 and 8 CP ) depend on eleven parameters: the three relevant oscillation parameters 
Am^, 612, di3 and the eight reduced solar fluxes /j. With the data from the different 
data samples (D) and the theoretical predictions for them in terms of these parameters 
uj = (Am^, 0i2, 613, /pp, • • • , /hep) we build the corresponding likelihood function 



£(D|w) = ^exp 



(2.2) 



where is a normalization factor. In Bayesian statistics our knowledge of uj is summarized 
by the posterior probability distribution function (p.d.f.) 

where Tr(uj\V) is the prior probability density for the parameters. In our model independent 
analysis we assume a uniform prior probability over which we impose the following set of 
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constraints to ensure consistency in the pp-chain and CNO-cycle, as well as some relations 
from nuclear physics: 

• The fluxes must be positive: 

$i>0 => fi>0. (2.4) 



• The number of nuclear reactions terminating the pp-chain should not exceed the 
number of nuclear reactions which initiate it [30,48]: 

$7 Bc + $8 B < $pp + $pcp 

=► 8.49 x 10- 2 Abc + 9.95 x KT 5 /«b < / PP + 2.36 x 10- 3 / P c P ■ (2.5) 

• The 14 N(p, 7) 15 reaction must be the slowest process in the main branch of the 
CNO-cycle [48]: 

$15 < $13 N => /wo < L 34/l3 N (2.6) 

and the CNO-II branch must be subdominant: 

$17 F < $15 => /l7 F < 37/l5Q • (2.7) 



• The ratio of the pep neutrino flux to the pp neutrino flux is fixed to high accuracy 
because they have the same nuclear matrix element. We have constrained this ratio 
to match the average of the BPS08(GS) and BPS08(AGS) values, with la Gaussian 
uncertainty given by the difference between the values in the two models 2 

f 

= 1.008 ±0.010. (2.8) 

/pp 

Following standard techniques we reconstruct the posterior p.d.f. in Eq. (2.3) using a Monte- 
Carlo algorithm; full details of our approach are given in Appendix B. 

The number of independent fluxes is reduced when imposing the so-called "luminosity 
constraint", i.e., the requirement that the sum of the thermal energy generation rates 
associated with each of the solar neutrino fluxes coincides with the solar luminosity [49] : 

i=l 

Here the constant on is the energy provided to the star by the nuclear fusion reactions 
associated with the i th neutrino flux; its numerical value is independent of details of the 
solar model to an accuracy of one part in 10 4 or better [30]. A detailed derivation of this 
equation and the numerical values of the coefficients «j , which we reproduce for convenience 

2 We have verified that assuming a flat distribution over the la uncertainty interval does not produce 
significant differences in the results of our analysis. 
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in Table 1, is presented in Ref. [30]. In terms of the reduced fluxes Eq. (2.9) can be written 
as: 



l = J>/i with fa = 



i=l 



L /[4tt(A.U.)2 



(2.10) 



where fa is the fractional contribution to the total solar luminosity of the nuclear reactions 
responsible for the production of the <3?- ef neutrino flux, and L@/[4ir (A.U.) 2 ] = 8.5272 x 
10 11 MeVcm~ 2 s _1 [30]. The analysis performed incorporating the priors in Eqs. (2.4-2.9) 
will be named "analysis with luminosity constraint" , V = L G , and for this case the prior 
probability distribution is: 



n(u'\L(7 



1 

N 




exp 



(W/pp- 1-008) ; 

2<7 2 



if Eqs. (2.4-2.7) and (2.9) are verified, 

otherwise, 

(2.11) 

where A is a normalization factor and a = 0.010. When only Eqs. (2.4-2.8) are imposed 
we will speak of "analysis without luminosity constraint", V = $ , so: 



1 

A 




exp 



(W/pp - 1 -008) ^ 
2<r 2 



if Eqs. (2.4-2.7) are verified, 
otherwise. 



(2.12) 



Let us notice that the conditions in Eqs. (2.4-2.7) and Eq. (2.9) are constraints on some 
linear combinations of the solar fluxes and they are model independent, i.e., they do not 
impose any prior bias favouring either of the SSM's. Furthermore we have chosen to center 
the condition (2.8) at the average of the BPS08(GS) and BPS08(AGS) values, with la 
Gaussian uncertainty given by the difference between the values in the two models, to 
avoid the introduction of a bias towards one of the models. In the next sections we will 
comment on how our results are affected when this prior is centered about the BPS08(GS) 
or the BPS08(AGS) prediction. 



3. Results 

Our results for the analysis with luminosity constraint are displayed in Fig. 1, where we 
show the marginalized one-dimensional probability distributions p(/j|D,L Q ) for the eight 
solar neutrino fluxes as well as the 90% and 99% CL two-dimensional allowed regions. 
The corresponding ranges at la (and at the 99% CL in square brackets) on the oscillation 
parameters are: 

Amli = 7.6 ± 0.2 [±0.5] x 10" 5 eV 2 , 
sin 2 0i2 = 0.33 ± 0.02 [±0.05] , (3.1) 
sin 2 6»i3 = 0.02 ± 0.012 , 
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Figure 1: Constraints from our global analysis on the solar neutrino fluxes. The curves in the 
rightmost panels show the marginalized one-dimensional probability distributions, before and after 
the inclusion of the Borexino spectral data. The rest of the panels show the 90% and 99% CL 
two-dimensional credibility regions (see text for details). 
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As mentioned above we have checked the stability of the results under changes in the 
assumption of the prior in Eq. (2.8). We find that if we center this prior at the BPS08(GS) 
prediction (f pep / f pp = 1) the best fit value for pep neutrinos is changed to f pep = 0.986 
(c&pep = 1.390 x 10 8 cm _2 s _1 ). Conversely if the Gaussian prior in Eq. (2.8) is centered 
at the BPS08(AGS) prediction (f pcp /f pp = 1.016) we get f pep = 1.003 ($ pep = 1.414 x 
10 8 cm _2 s _1 ). All other fluxes are unaffected. 

For the sake of illustration we have also performed a Gaussian fit to the two-dimensional 
p.d.f. for the eight fluxes. The best Gaussian approximation to the real p.d.f. is charac- 
terized by flux uncertainties obtained by symmetrizing the la ranges quoted in Eq. (3.2), 
and by the following error correlation matrix: 
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(3.3) 



As seen in Fig. 1 and in Eq. (3.3) the most important correlation appears between the pp 
and pep fluxes, as expected from the relation (2.8). The correlation between the pp (and 
pep) and 7 Be flux is directly dictated by the luminosity constraint (see comparison with 
Fig. 3). All these results imply the following share of the energy production between the 
pp-chain and the CNO-cycle 

£pp-chain no ^+0.005 [+0.011] , , ^CNO _ n nl , +0,006 [+0.014] (o A \ 

_ U.986_ 006 [_ aol4 J — = - 0.014_ 005 Uo.OllJ ' K 6A ) 

in perfect agreement with the SSM's which predict Lcno/^o < 1% at the 3a level. 

The sensitivity of the various experiments is illustrated in Fig. 2, where we plot the 
contribution of each flux to the total event rates at the radiochemical experiments as well as 
SNO and SK (for Borexino see Appendix A) together with the corresponding experimental 
values and uncertainties. To highlight the sensitivity to the hep flux we plot separately 
the rate for the last energy bin in SK (SK-hi, E e > 16 MeV); similar results hold for 
the highest energy bins of SNO. The rates are computed for the best fit value of the 
oscillations parameters, Eq. (3.1). We show the contributions as predicted by BPS08(GS) 
and BPS08(AGS) solar models and by our best fit values for the fluxes given in Eq. (3.2). 

In order to check the consistency of our results we have performed the same analy- 
sis without imposing the luminosity constraint, Eq. (2.9). The corresponding results for 
p(fi\D, ]Jjq) and the two-dimensional allowed regions are shown in Fig. 3. As expected, the 
pp flux is the most affected by the release of this constraint. This is so because the pp 
reaction gives the largest contribution to the solar energy production, as can be seen in 
Table 1. Imposing the luminosity constraint as an upper bound on the pp flux would imply 
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Figure 2: Contribution of each solar neutrino flux to the total event rates at different experiments. 
The oscillation parameters are set to their best fit value, Eq. (3.f). We show the contributions as 
predicted by BPS08(GS) and BPS08(AGS) solar models as well as our best fit values given in 
Eq. (3.2). 



that this flux cannot exceed its SSM prediction by more than 9%. Conversely, releasing 
this constraint allows for a much larger pp flux. The pep flux is also severely affected due 
to its strong correlation with the pp flux, Eq. (2.8). On a smaller scale the CNO fluxes 
are also affected, mainly as an indirect effect due to the modified contribution of the pp 
and pep fluxes to the Gallium and Chlorine experiments, which leads to a change in the 
allowed contribution of the CNO fluxes to these experiments. Thus in this case we get: 

f _ f _ n Qs+ - 16 r+0-47i 

Jpp — Jpep — U.»O_ 15 [_ .40J ' 

f, _ i ni+ - 1 [+°- 27 i 

J?Be — 1 - ul -0.09 L-0.22J ' 

/xsn = 2.7±1J [±H] , (3.5) 

/i5 = 1.9 ± 1 [±?:i] , 

/ir F < 34 [79] . 

The determination of the 8 B and hep fluxes (as well as the oscillation parameters) is 
basically unaffected by the luminosity constraint. 

Interestingly, the idea that the Sun shines because of nuclear fusion reactions can 
be tested accurately by comparing the observed photon luminosity of the Sun with the 
luminosity inferred from measurements of solar neutrino fluxes. We find that the energy 
production in the pp-chain and the CNO-cycle without imposing the luminosity constraint 
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Figure 3: Same as Fig. 1 but without the luminosity constraint, Eq. (2.9). The curves in the 
rightmost panels show the marginalized one-dimensional probability distributions, before and after 
the inclusion of the Borexino spectral data. The rest of the panels show the 90% and 99% CL 
two-dimensional credibility regions (see text for details). 



are given by: 
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(3.6) 



Comparing Eqs. (3.4) and (3.6) we see that the luminosity constraint has only a limited 
impact on the amount of energy produced in the CNO-cycle. However, as discussed above, 
the amount of energy in the pp-chain can now significantly exceed the total quantity allowed 
by the luminosity constraint. Altogether we find that the present value for the ratio of the 



10 



neutrino-inferred solar luminosity, Lq (neutrino-inferred), to the photon luminosity L & is: 

= 1.00 ±0.14 [+° f 4 ] . (3.7) 



L (neutrino-inferred) = ^ ± QU ^ 



Thus we find that, at present, the neutrino-inferred luminosity perfectly agrees with the 
measured one, and this agreement is known with a la uncertainty of 14%. 

3.1 The role and potential of Borexino 

As seen in Figs. 1 and 3 the inclusion of Borexino has a very important impact on the 
determination of the 7 Be, pep and CNO fluxes, and indirectly on the pp flux. As mentioned 
above and described in Appendix A, in our analysis we have fitted the 160 data points of the 
Borexino energy spectrum in the 365-2000 keV energy range, leaving as free parameters 
the normalization of the the 11 C, 14 C, 210 Bi and 85 Kr backgrounds, the three relevant 
oscillation parameters, and the normalization of all the solar neutrino fluxes. In contrast, 
the Borexino collaboration fits the spectrum in the full energy range 160-2000 keV, and 
allows for free normalizations of the n C, 10 C, 210 Bi+CNO and 85 Kr backgrounds as well as 
14 C (which introduces an overwhelming background but is only relevant for events below 
250 keV, hence it does not contribute to our analysis). Besides the normalization of these 
background components only the 7 Be flux normalization is fitted to the data, and no direct 
information on the normalization of the other solar fluxes is extracted. In particular, the 
CNO fluxes are added to the 210 Bi background and fitted as a unique "background" , while 
the other solar fluxes are fixed to the BPS08(GS) prediction and the oscillation parameters 
are fixed to the best fit point of the global pre-Borexino analysis. With this procedure 
they determine the interaction rate for the 0.862 MeV 7 Be neutrinos to be 49±3 s t a t ±4 sys , 
which corresponds to fv^ e = 1-02 ± 0.10. Given the precision of the data and the energy 
spectrum of the irreducible backgrounds, their procedure is perfectly acceptable for the 
purpose of extracting the 7 Be normalization alone. However, in this work we are interested 
in testing the full set of SSM fluxes, not just 7 Be. In this case the consistent procedure is 
to allow for independent normalizations of all the solar fluxes. 

In order to illustrate the impact of these two different approaches on the determination 
of the SSM fluxes, we have repeated our analysis using as unique data input from Borexino- 
LE the total interaction rate for 7 Be neutrinos quoted by the collaboration. The results 
are shown in Fig. 4. As seen in this figure the constraints imposed on the 7 Be flux after 
the inclusion of Borexino are basically the same, irrespective of whether one includes the 
full Borexino energy spectrum in the range 365-2000 keV or just the total 7 Be event rate 
extracted by the Borexino collaboration. However, the best fit and allowed ranges for 
the CNO fluxes are not the same. This proves that, despite the unknown level of 210 Bi 
contamination, the Borexino spectral data can still provide useful information on the CNO 
fluxes. When using only the total 7 Be event rate this information is lost and the constraints 
on CNO arise exclusively from the Gallium and Chlorine experiments. 

As shown in Fig. 4, the inclusion of the complete Borexino-LE spectrum leads to an 
improvement (albeit not very significant) of the determination of the 13 N flux. Without 
Borexino this flux is mostly (and poorly) constrained by the Gallium experiment, and 



- 11 - 



O 




0.9 0.95 1 1.05 




0.5 0.75 



1.25 1.5 



O 

n 





O 



CL 20 




0.9 0.95 1 1.05 




0.5 0.75 1 1.25 1.5 




J3 1 
O 

CL 0.5 



: M I 1 I I 


i i i i i 
T 150 = 








i i i i i 







4 6 8 



0.5 



I J I I I I I I I 

— - 


I I I J I I ! 


t 




^^^i^^sw I I 






i i i i i 



8 10 



Figure 4: Marginalized one-dimensional probability distributions for the fluxes contributing to 
Borexino-LE. The full line shows the determination obtained by fitting the Borexino spectrum data 
as described in Appendix A. The dashed-dotted line shows what the results would be if instead one 
had used the Borexino result for the extracted interaction rate of the 0.862 MeV 7 Be neutrinos. The 
dotted line represents the precision obtainable with the simulated "ideal" spectrum as described in 
Sec. 3.1. 
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including the additional information from Borexino positively adds to its knowledge. We 
notice, however, that the best fit value of the 13 N flux in either analysis is always higher 
than the prediction of any of the SSM's (although fully compatible at better than 1.5a). 
This behavior is driven by the Gallium rate which is slightly higher than expected in any 
of the SSM's within the framework of three neutrino oscillations. A higher best fit value 
of 13 N can easily accommodate this observation without conflicting with any of the other 
experiments nor with the observed spectrum at Borexino. On the contrary, this is not 
the case for the 15 flux. Adding the information from Borexino-LE spectrum leads to a 
(also small) worsening of the precision in the determination of this flux. We traced this 
apparently counter-intuitive result to the existing tension between the low Chlorine rate and 
the predicted rate within the framework of three neutrino oscillations. As a consequence of 
this tension, Chlorine pushes towards lower values of the 15 flux, whereas the spectrum 
of Borexino-LE prefers a higher amount of 15 0. When only the total event rate of Be is 
used in the fit the extracted 15 O flux is mostly driven by the Chlorine result. When the 
Borexino-LE spectrum is also included the tension results into a higher best fit for the 15 
flux and a worsening of the precision. 

Among all solar neutrino fluxes, the CNO ones are those for which the differences 
between the BPS08(GS) and BPS08(AGS) SSM predictions are largest. It is therefore 
interesting to explore whether this discrepancy can be resolved with future Borexino data, 
and to what degree the CNO and pep fluxes can be better determined. Besides the accu- 
mulation of more statistics, in the near future Borexino aims at reducing the systematic 
uncertainties with the deployment of calibration sources in the detector [40]. Ideally, if 
the 11 C background could be subtracted from the signal the pep and CNO fluxes would 
become directly accessible. In order to illustrate the potential of this perspective we have 
simulated an "ideal" spectrum of 85 bins in the energy range 365-1238 keV according to 
the expectations from the central values of the BPS08(GS) fluxes and the best fit point of 
oscillations. In our simulation we have assumed that the n C has been fully removed, while 
the other backgrounds are added under the same assumptions as in the present Borexino 
analysis. We have also assumed double statistics and a reduction by a factor three of the 
systematic uncertainties. The results of this fit are presented in Fig. 4. As can be seen, 
with this ideal experiment the level of accuracy can be substantially improved for most 
fluxes, with the exception of 13 N flux whose determination becomes less precise. This is 
a consequence of the tension between the higher value of 13 N preferred by the Gallium 
experiments and the SSM value used for the simulated spectrum: if we had simulated data 
corresponding to a higher value of 13 N flux, the precision in the determination of this flux 
would also have improved. In any case, our results show that even with this optimistic 
improvement of Borexino the precision of the CNO fluxes remains far below the present 
uncertainties of the SSM's. 

3.2 Comparison with the Standard Solar Model(s) 

In Fig. 5 we show the marginalized one-dimensional probability distributions for the solar 
neutrino fluxes as determined by our analysis, together with the predictions for the two 
SSM's given in Ref. [9]. In order to statistically compare our results with the SSM's 



- 13 - 




O 



4 |T 

3.2 V 
2.4 f 
1.6 h 
0.8 I 






150 



1 



15 
12 
9 
6 
3 



Jill 


i i i i i 1 i i i 1 i_ 

> Fit = 

BPS08(GS)z 

BPS08(AGSJ 


Z / 1 

/ 

_ / 





0.6 0.8 1 1.2 1.4 



Figure 5: Marginalized one-dimensional probability distributions for the best determined solar 
fluxes in our analysis as compared to the predictions for the two SSM's in Ref. [9]. 



predictions we perform a significance test. We start by constructing a posterior probability 
distribution function for the solar fluxes as well as for the SSM central values. This is 
defined as the probability distribution from the data subject to the prior distribution of an 
arbitrary SSM: 

„(/, £-|D, SSM) _ ^l/W./^'ISSM) 



where 



-2 In 



fC(B\f) vr(/, /SSMJSSM) d/d/SSM 



n(f, /^ISSM)] = J2(fi ~ /7 SM )^ssm,^ " /T SM ) (3-9) 

hi 



and Vssm is the covariance matrix for the assumed SSM model. We build the covariance 
matrix for arbitrary models by interpolating the covariance matrices for the BPS08(GS) 
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Figure 6: Probability distribution function g(t) (full line, see text for details). For comparison we 
show the corresponding distribution for a \ 2 P-d.f. with 8 degrees of freedom (dashed line). 



and BPS08(AGS) models given in Ref. [9]. 3 Since these covariance matrices are very similar 
to each other, our results are not sensitive to this assumption. Furthermore, in order to 
improve this approximation one would need a continuous model dependence for the flux 
covariance matrix, which is currently unavailable. 

The posterior probability distribution for a SSM characterized by given central values 
and covariance matrix, subject to the constraints imposed by the data is then 

p(f SM \B) = J p(f,f SM \D,SSM)df (3.10) 

Fromp(/ SSM |D) we define a probability distribution function for the statistics t as 

g(t) = J P (F SM \D)5[t + 2\n(p(^ SM \D))] df^ SM (3.11) 

By definition g(t) is a function normalized to 1 in the interval t m j n < t < oo. With this 
definition, t would follow a \ 2 distribution if p(/|D,Lq) were exactly Gaussian. In Fig. 6 
we plot the function g{t). For comparison we show the corresponding distribution for a \ 2 
p.d.f. with 8 degrees of freedom (dashed line). 

The significance of the agreement between the data and what is expected under the 
assumption of a given model is quantified in terms of the probability P agr , defined as the 



3 The correlation matrix is given in http://www.mpa-garching.mpg.de/~aldos. 



- 15 - 



probability to find t in the region of equal or larger compatibility with the data than the 
level of compatibility observed within the given model: 

r ^max 

^g S s(ags)=/ 9{t)dt (3.12) 

^*GS(AGS) 

where tcs(AGS) is the value of the statistic obtained for the central value fluxes of the 
specific model: 

^GS(AGS) = -2 In [p(/ GS(AGS) |D)] . (3.13) 

We found that the GS model has a lower t, t G s = 8.5, while £ags = 110. With the 
probability distribution g(t) shown in Fig. 6 this corresponds to = 43% and -P^GS = 
20%. 

For comparison we have also constructed a \ 2 function comparing the best fit values 
of the fluxes in each of the models with those obtained in the analysis without the SSM 
priors and with the uncertainties given by the combined covariance matrix 



X 2 = £Uf S(AGS) - IP) fe(AGS) + (if (AGS) - I?) • (3.14) 

ij 



Here Vb is the covariance matrix obtained by the best Gaussian approximation to the 
p(fi\D, Lq) probability distribution function, Eq. (3.3), and ff are the best fit fluxes from 
the data analysis without any SSM prior, Eq. (3.2). If the distribution p(/|D,L ) were 
exactly Gaussian, both tests would be equivalent. We found that this test still yields a 

better fit for the GS model, x GS = 5 - 2 ( P GS = 74% ) and Xags = 5 - 7 ( P ags = 68% )> 
but gives a higher probability for both models. We also find that if the prior in Eq. (2.8) 
is centered at the BPS08(AGS) prediction the analysis renders the same probability for 
both models. Conversely if Eq. (2.8) is centered at the BPS08(GS) prediction the slight 
preference for the BPS08(GS) model found above is enhanced to Xqs = ^-2 (P«l r = 74%) 
versus X | GS = 7.4 (P*g s = 50%). 

From these results we conclude that, while the fit shows a slightly better agreement 
with the GS model corresponding to higher metallicities, the difference between the two 
is not statistically significant. This is partly due to the lack of precision of present data. 
But we also notice that, while the measurements of SNO and SK favor a lower 8 B flux as 
predicted by the low metallicity models, the determination of the 7 Be flux in Borexino and 
the corresponding determination of the pp flux from the luminosity constraint show better 
agreement with the GS predictions. 



4. Summary 

We have performed a solar model independent analysis of the solar and terrestrial neu- 
trino data in the framework of three-neutrino oscillations, following a Bayesian approach 
in terms of a Markov Chain Monte Carlo using the Metropolis-Hastings algorithm. This 
approach has allowed us to reconstruct the probability distribution function in the entire 
eleven-dimensional parameter space, consistently incorporating the required set of theoret- 
ical priors. The best fit values and allowed ranges for the eight solar neutrino fluxes are 
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summarized in Eq. (3.2) and Fig. 1 for the analysis with the luminosity constraint, and in 
Eq. (3.5) and Fig. 3 for the more general case of unconstrained solar luminosity. We found 
that at present the neutrino-inferred luminosity perfectly agrees with the measured one 
and it is known with a la uncertainty of 14%. We have also tested the fractional energy 
production in the pp-chain and the CNO-cycle, finding that the total amount of the solar 
luminosity produced in the CNO-cycle is bounded to be Lcno/Ad < 2.8% at 99% CL 
irrespective of whether the luminosity constraint is imposed or not. 

We have then presented a statistical test which can be performed with these results in 
order to shed some light on the so-called solar composition problem, which at present arises 
in the construction of the Standard Solar Model. We found that the low value of the 8 B flux 
measured at SK and SNO points towards low metallicity models, whereas the measurement 
of Be in Borexino and the corresponding value of the pp flux implied by the luminosity 
constraint show better agreement with high metallicity models. Altogether the fit shows a 
slight preference for models with higher metallicities, however the difference between the 
two models is not very significant at statistical level. While a realistic improvement of the 
Borexino data analysis in the near future can positively affect the direct determination of 
most solar neutrino fluxes, it is unlikely that enough precision will be achieved to go beyond 
the the present theoretical uncertainties of the SSM's. The largest difference between the 
models lies on the CNO fluxes that give predictions which differ by about 30%. Thus 
ideally in order to achieve a statistically meaningful discrimination between the models 
one would need a low energy solar neutrino experiment capable of measuring the neutrino 
energy spectrum for energies between 0.5 MeV < E v < 1.5 MeV and, more importantly, of 
rejecting the radioactive backgrounds to the required level, so to allow for a determination 
of the CNO fluxes at 0(30%) precision. 
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A. Analysis of Borexino spectra 

In our analysis of Borexino we include both the low-energy (LE) data presented in Ref. [40], 
which are crucial for the reconstruction of the 7 Be line, as well as the high-energy (HE) data 
discussed in Ref. [38], which are mostly sensitive to the 8 B flux. For the low-energy part we 
extracted the 180 experimental data points and the corresponding statistical uncertainties 
from Fig. 2 of Ref. [40], checking explicitly that the statistical error o| tat is just the square 
root of the number of events N? x in each bin b (except in the region where statistical a's 
subtraction had been performed). Similarly, for the high-energy part we extracted the 6 
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experimental data points and statistical uncertainties from Fig. 3 of Ref. [38]. For both 
data sets the theoretical prediction A^ h for the bin b is calculated as follows: 

~dK iK ^ w e {Eu > Te) Rb{Te & dEu + (A - 1} 

where uj describes both the neutrino oscillation parameters and the eight solar flux normal- 
izations, and |* is a set of variables parametrizing the systematic uncertainties as required 
by the "pulls" approach to % 2 calculation. Here n e \ is the number of electron targets in a 
fiducial mass of 78.5 tons with an electron/nucleon ratio of 11/20 for pseudocumene, and 
T b run is the total data-taking time which we set to 192 and 246 live days for LE and HE 
data, respectively. In the previous formula do~ a /dT e is the elastic scattering differential 
cross-section for neutrinos of type a £ {e, fi, r}, and d$^ et / dE v is the corresponding flux of 
solar neutrinos at the detector - hence it incorporates the neutrino oscillation probabilities. 
The detector response function Rb(T e \£) depends on the true electron kinetic energy T e and 
on the three systematic variables £ vo i, £ S ci and £ res : 

/•T-ax (1+7r 6 ci £ scl) 

R b (T e \0 = (1 + vr vol £voi) / Gauss [T e - T , a T (l + 7r res &«,)] dT' (A.2) 

where Gauss(x, a) = exp [— x 2 /2a 2 ] j\pli:a is the normal distribution function, while T h mm 
and T 6 max are the boundaries of the reconstructed electron kinetic energy T' in the bin 
b. Note that we assumed an energy resolution ax/T e = 6%/y / T e [MeV], rather than the 
"official" value 5%/ \/T e [MeV] quoted by the collaboration, since our choice lead to a 
perfect match of the Be line shown in Fig. 2 of Ref. [40] . We verified that also the other 
solar fluxes plotted in Ref. [50] are carefully reproduced. As for the effects introduced by 
systematic uncertainties, we assumed 7r vo i = 6% for the fiducial mass ratio uncertainty, 
7Tg cl = 2.4% (1%) for the energy scale uncertainty in LE (HE) data, and an arbitrary 
7r res = 10% for the energy resolution uncertainty. 

The backgrounds A^ kg (£) which appear in Eq. (A.l) only affect the low-energy data, 
and are not included in the calculation of the high-energy event rates. The 10 C, 11 C, 14 C 
and 85 Kr background shapes were taken from Fig. 2 of Ref. [40], whereas the 238 U, 214 Pb 
and 210 Bi were extracted from slide 7 of Ref. [50]. We explicitly verified that with the 
normalizations as inferred from these figures the sum of all these backgrounds with the 
expected SSM fluxes precisely reproduces the "Fit" line shown in Fig. 2 of Ref. [40]. Note 
that, because of the overwhelming 14 C at energies below ~ 250 keV, our fit in this energy 
region is never good. We do not know if this is due to the loss of numerical precision in 
our extraction of the 14 C shape or to the fact that there are additional free parameters to 
be fitted for this background. In any case, in order to avoid biasing our analysis by the 
low quality description on these data points we use only the 160 points of the spectrum 
above 365 keV. Hence the 14 C background is irrelevant. Following the procedure outlined 
in Ref. [40] the normalization of the 238 U and 214 Pb backgrounds are assumed to be known, 
whereas the normalizations of the 85 Kr, 210 Bi, 11 C and 10 C backgrounds are introduced as 
free parameters and fitted against the data - taking care to ensure their positivity. Hence: 

» r bkg/A _ A/-U238 , ArPb214 , i\rKr85 f , ArBi210/; , ivrCllt , AT CW C (\ 1\ 

(?) — + iV f> + 7V 6 ?Kr85 + ^b ?Bi210 + ^b SiCll + iV fe ?C10 • (A. 3) 
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Figure 7: Spectrum for the best fit point of our spectral fit to the Borexino-LE data in the energy 
region between 350-2000 KcV under the assumptions described in the Appendix (left) and A% 2 as 
a function of the 7 Be flux for the different test analysis of Borexino data (right). 



The x 2 (3) function is constructed in the usual way in the context of the pull method, 
by introducing standard penalties for the £ variables (except for those parametrizing the 
free normalizations of the backgrounds) and marginalizing over them: 



X 2 (w) = min 



,,.stat 



+ £vol + 4 2 cl+£res> • (A.4) 



As a test of our procedure we first perform a fit under the same assumptions as Ref. [40], 
i.e., besides the backgrounds we only fit the Be flux normalization to the data. The other 
solar fluxes are fixed to their BPS08(GS) prediction and the oscillation parameters are fixed 
to the best fit point of the global pre-Borexino analysis. The results of this test are shown in 
Fig. 7. Comparing the left panel with Fig. 2 of Ref. [40] we observe a perfect agreement in 
the best fit 7 Be flux spectra. In the right panel we plot the A% 2 for this test fit as a function 
of /7Be) f° r both Borexino-LE alone and the combination of Borexino-LE and Borexino- HE 
data. As can be seen, for Borexino-LE our procedure leads to a determination of the 7 Be 
normalization in very good agreement with the value JVgg = 1.02 ± 0.10 obtained by the 
Borexino collaboration. The inclusion of the Borexino-HE tends to push the extracted 
value of /7 Be towards a slightly higher value. This is due to the assumed correlation of the 
systematic uncertainties (in particular the one associated with the total fiducial volume) 
between LE and the HE. This small effect is diluted once both data sets are included in 
the global fit, and the final results are practically independent of the degree of correlation 
assumed between the systematic errors. 



B. Details of the Markov Chain Monte Carlo 

In this analysis we have used the Metropolis-Hasting algorithm including an adapting 
algorithm for the kernel function to increase the efficiency. The algorithm is defined as 
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follows: 

1. Given a parameter set Co, a new value Co' is generated according to a transition kernel 
q(Co,Co'). We start with a flat kernel and, once the chain has reached a certain size, 
we use a kernel in terms of the covariance matrix V computed with the points in the 
chain. If U is the matrix diagonalizing V and di are the eigenvalues, Co' = Co + U Co 
with Cji generated according to a distribution x exp(— Cbi/di). The kernel is 
adapted, i.e., the covariance matrix is recalculated, every several steps. 

2. With 63 and uj' we compute the value: 

3. A random number < r < 1 is generated and if r < h, Co' is accepted in the chain. 

4. We go back to step (1), starting with Co' if it has been accepted or again with Co if 
not. 

All the points accepted in this algorithm constitute the Markov Monte-Carlo Chain {Co a } 
with a = 1, . . . ,N to t, where N tot is the total number of points in the chain. The method 
ensures that, once convergence has been reached, the chain takes values over the parameter 
space with frequency proportional to the posterior probability distribution function. 

Technically, in order to reconstruct the posterior p.d.f. from the chain we discretize 
the parameter space by dividing the physically relevant range of each parameter uj{ into 
rii subdivisions of length £^ (with 1 < ki < rij). Denoting as £lki...k m the cen corre- 
sponding to subdivisions A^ 1 . . . A^ 71 (m = 10 or 11 for the analysis with or without the 
luminosity constraint, respectively), we compute the value of the posterior p.d.f. as 

p(Co G n kl ... km \D,V) = %^ (B.2) 

where M^...^ is the number of points in the chain with parameter values within the cell 
^fci...fc m > and Vfci...fc m = ^\ x • • • x l^ 1 is the volume of the cell. In order to ensure that the 
procedure generates a smooth p.d.f., a sufficiently large N tot is needed. 

The marginalized one-dimensional p.d.f. for the parameter ooi is reconstructed as 

Tlj 

p(^GAMD,P) 1 . dim = i- £ A/ y-*" ■ (B.3) 

H k^i=i tot 

Similarly, the marginalized two-dimensional p.d.f.'s for the parameters (ooi,ooj) is 

P ( Wi €A? >Wi €A>,P) Mtal = ' E M ^-^ , (B.4) 

From these, we obtain the two-dimensional credibility regions with a given CL as the region 
with smallest area and with CL integral posterior probability. In practice they are obtained 
as the regions surrounded by a two-dimensional isoprobability contour which contains the 
point of highest posterior probability and within which the integral posterior probability 
is CL. 
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